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ABSTRACT 

Vertical  Roller  Mill  (VRM)  has  been  promptly  employed  and  extensively  used  for  clinker  and  cement  grinding 
for  past  35  years,  because  ofits  higher  efficiency  and  low  power  consumption.  VRM  is  a  highly  nonlinear,  time  delay  and 
composite  industrial  process  with  several  process  loops  having  a  robust  coupling  between  them.  Time  series  data  is 
collected  from  a  real  time  cement  plant  with  grinding  capacity  of  170-300  tph  for  system  identijication  purpose. 
The  collected  industrial  data  is  normalized  and  preprocessed  using  a  moving  average  window  filter  to  remove  odd 
samples;  trends  and  means  are  removed  for  outlier  elimination  and  missing  values,  then  the  data  is  divided  into  an 
estimation  data  and  a  validation  data.  In  this  paper  linear  and  nonlinear  system  identification  methods  are  applied  to 
identify  the  VRM  model.  Tinally,  this  paper  proposes  an  accurate  VRM  model  for  industrial  requirement  to  adopt  modern 
control  techniques  like  Model  Predictive  Control  (MPC)  and  Adaptive  Control  (AC)  for  cement  grinding  in  order  to 
improve  its  automatic  control  level. 

KEYWORDS:  Vertical  Roller  Mill,  System  Identification,  Hammerstein  Model,  Polynomial  Model  &  Auto  Regressive 
Models 
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1.  INTRODU CTION 

VRMs  have  evolved  through  the  years  and  are  currently  being  used  at  various  locations  worldwide  for 
grinding  slag,  OPC  grade  43  and  53,  PPC  and  PSC  cement.  Clinker  and  cement  grinding  machines  are  the 
maximum  energy  consumption  units  in  the  cement  manufacturing  process.  Figure  1.  Shows  the  various  parts  of  the 
VRM  used  for  cement  grinding.  Figure  2.  Shows  the  snapshot  of  the  real  time  VRM  used  for  cement  grinding.  The 
working  principle  of  VRM  involves  start,  steady  state  stop  processes.  The  important  processes  present  in  VRM  are 
grinding,  drying  and  separation. 
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Figure  1:  Parts  of  Vertical  Roller  Mill  Adopted  from 
https://www.scribd.com/doc/97749158/Vrm-Presentation 


A  collective  appliance  of  pressure  frame,  hydraulic  tensioning  system  and  pull  rod  is  used  to  force  the  rollers 
against  the  grinding  bed.  From  the  air  stream  the  completed  product  is  detached  and  conveyed  to  the  cement  silos  [1],  [2] 
and  [6]. 


Figure  2:  Photo  of  the  Real  Time  VRM 
Used  in  Cement  Plant 


2.  MODEL  IDENTIFICATION 

The  model  identification  is  the  development  and  investigation  of  the  models  in  accordance  with  the  input-output 
data.  Models  are  built  from  basic  physical  laws  and  other  well  recognized  relationships.  The  flow  chart  of  the  system 
identification  loop  using  input-output  data  is  shown  in  Figure  3.  The  extensively  used  linear  models  are  transfer  function 
models,  Linear  Auto  regressive  (LARX)  models,  State  space  models  and  polynomial  models.  The  time  series  input-output 
data  is  the  basis  for  system  identification  irrespective  of  the  model  structure  [3].  The  generalized  model  structure  and  its 
predictor  are  represented  and  expressed  as  (1)  and  (2). 

A(p)*y(t)  =  p~"ky(t>  ^^*u(t)+^pl*e(t)  (1) 

F(p)  D(p) 
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y(t\0)  =  D(P)B(P)  u(t)  + 
C(p)F(p) 


1- 


D(p)A(p) 

C(p) 


y(t) 
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(2) 


Where,  u(t)~  cause,  y(t)~  response,  e(t)~  disturbance,  A(p),  B(p),  C(p),  D(p),  and  F(p)~  matrices  with  unknown 
parameters  p. 


Figure  3:  Flow  Chart  of  System  Identification 

3.  LINEAR  ESTIMATION  METHODS 

The  mathematical  model  commonly  expresses  the  dependencies  of  yarious  process  variables.  Transfer  function 
model  is  the  most  widely  used  methods  for  modelling  dynamic  linear  time-invariant  systems.  The  transfer  function  models 
are  discussed  in  [4],  [5]  and  [7]. 

3.1.  Transfer  Function  Models 

Transfer  function  based  mathematical  models  are  more  appropriate  for  the  time  domain,  frequency-domain  and 
stability  analysis  of  the  systems.  Continuous  and  a  discrete  time  systems  are  expressed  as  (3)  and  (4). 


y(t)  =  G(s )  *  u(t)  +  e(t) 
y(k)  =  G(z)*u(k)  +  e(k) 

Where,  G(s)  and  G(z)~  transfer  function  between  the  cause  and  the  effect. 


(3) 

(4) 


The  transfer  function  models  for  single  input-  single  output  continuous  time  and  discrete  time  systems  are 
expresses  as  (5)  and  (6). 


G(s)  = 


G(z)  - 


fr0+V  +  ---+£m_1.ym'  +  frms 
a0  +aj5H — an_vsn~]  +ansn 


ao  +alz  +  ---an_lzn~l  +anzn 


(5) 


(6) 


Where,  bm  and  an-  coefficients  of  the  numerator  polynomial  function,  m  and  n-  order  of  the  numerator  and 
denominator. 


The  equations  for  the  continuous  and  discrete  transfer  function  of  multi  input  single  output  systems  are 
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n 


yi  =HGij^ui 


(7) 


n 


yt  =HGij(^uj 


(8) 


Where,  Gir  transfer  functions  between  the  cause  and  the  effect,  I  and  j-  number  of  inputs  and  outputs. 

The  algorithmic  steps  to  construct  transfer  function  model  is  as  follows 
Stepl:  Apply  Prediction  Error  Method  to  calculate  discrete  model 
Step2:  Convert  discrete  model  to  a  continuous  model  by  applying  Zero-Order-Hold 
Step3:  Optimize  continuous  model  by  using  Gauss-Newton  method 
3.2.  Polynomial  Models 

Polynomial  models  separately  define  the  stochastic  and  deterministic  parts  of  a  system  where  as  transfer  function 
models  define  only  the  deterministic  part.  The  polynomial  models  are  discussed  in  [8].  Different  polynomial  models 
widely  used  for  system  identification  include  ARX  models,  Auto  regressiye  Moying  Ayerage  models,  Instrument  yariable 


models,  Box-Jenkins  Model  and  Output  error  model. 


3.2.1.  ARX  Modelling 


Figure  4  shows  the  ARX  model  block  diagram.  The  linear  difference  equation  represents  the  input-output 


relationship  as 


y(0  +  at  *  vV-l)  + . +a„o  *y(t-na)  = 


(9) 


b{  *u(t- 1)  + . +bnb  *u(t -nb)  +  e(t) 


Where,  0  =  [at  a2  ■■■  an  bl---bn  ]T  is  the  adjustable  parameter, 


(10) 


e 


I 


1/A 


B/A 


Figure  4:  The  ARX  Model  Structure 
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3.2.2.  Instrument  Variable  (IV)  Modelling 

In  this  modeling,  yariance  and  Mean  Square  Error  (MSE)  depend  on  the  choice  of  instruments.  If  the  noise  and 
instruments  are  correlated  or  model  order  selection  is  wrong  or  filter  dynamics  cancels  plant  dynamics  then  consistency  is 
lost.  The  IV  model  is  shown  in  Figure  5. 


e 


3.2.3.  ARMAX  Modelling 

In  ARX  model  the  properties  of  the  noise  term  is  not  more  adequate.  In  ARMAX  model,  moving  average  of  white 
noise  relates  the  noise  model  for  more  flexibility.  The  model  equations  are 

y(t)  +  al*y(t-\)  + . +  aHa  *y(t-na)  = 

bl*u(t-l)  + . +  bHb  * u(t - nb)  +  e(t)  +  (11) 

cx*e(t-  1)H - \-cn  *e(t-nc) 

By  assuming  C(p)  =  1  +  C{p~l  H - \-  Cncp~nc 


Equation  (11)  is  written  as 

A(p)*y(t)=B(p)*u(t)+C(p)*e(t)  and  corresponds  to 


G(P,d) 


B(p) 

A(p) 


and  H(p ,  G )  =  — 

A(p) 


(12). 


The  ARMAX  model  structure  is  shown  in  Figure  6 


e 


Figure  6:  The  Equation  Error  Model 
Family  Model  Structure 


3.2.4.  Output  Error  (OE)  Modelling 

In  OE  model  structure  the  G  and  H  transfer  function  denominators  have  the  term  A  as  a  mutual  factor.  Figure  7 
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shows  the  structure  of  OE  model.  Let  the  undisturbed  output  w,  then  the  model  equations  are  as  follows 

W)+fi  *Mt-nf)  = 

b{  *u(t— 1)4 - \-bUb  *u(t-nb) 

y(t)  =  w(t)  +  e(t) 


For,  F(q)  =  1  +  fq  1  H - 1-  fn  q  ,  the  model  and  the  parameter  vectors  are 


y(t)  =  ^^-*u(t)  +  e(t) 
F(p) 


bl  b2  •"  K  /l  /•••  fn, 


=0 


(13) 

(14) 


(15) 


(16) 


u 


Fl)—  * 


Figure  7:  The  Output  Error  Model  Block  Diagram 
3.2.5.  Box-Jenkins  Model  Structure 


By  describing  the  OE  model  as  ARMA  model  gives 


y(t)  = 


^mo+^ 

F(p)  D(p) 


*e(t) 


(17) 


The  G  and  H  transfer  functions 


are  parameterized  as  a  rational  functions  independently.  Figure  8  shows  the  NJ 


model  structure. 


e 


Figure  8:  The  BJ  Model  Structure 

3.3.  State  Space  Modelling 

The  relationship  between  the  cause,  response  and  noise  is  written  as  difference  or  differential  equation  in  state 
space  modelling  by  using  an  auxiliary  state  vector  x(t).  The  model  equation  is  as  follows 
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x(t)  =  A(0)*  x(0  +  B{0)*u{t)  (18) 

Where,  A  and  B-  nxn  and  nxm  matrices. 

Let  m(t)  be  the  ideal  measurements  obtained,  m(t)=  H*x(t)  using  q  for  the  differentiation  operator  the  model 
equations  is  as  follows, 

[q*I-M0)]**(f)=B(0)*u(f)  (19) 

The  transfer  operator  from  u  to  m  is  m(t)  =  Gc  ( q ,  6)  *  u(t)  (20) 


Gc(q,0)  =  H*[q*I-  A(0)\l  * B(0) 

4.  NONLINEAR  MODELLING 


(21) 


A  non-linear  relationship  between  the  input  and  output  sequences  evidently  provides  superior  possibilities  to 
describe  a  system.  The  different  input  and  output  nonlinearities  are  wavelet  network,  piecewise  linear  function,  saturation, 
sigmoid  network,  one-dimensional  polynomial,  dead  zone,  or  a  custom  network.  Nonlinear  models  are  mainly  of  Nonlinear 
ARX  models  and  Hammerstein-Wiener  models. 


4.1.  Nonlinear  ARX  (NARX)  Modelling 

NARX  model  comprises  of  nonlinearity  estimator  and  model  regressors.  Both  linear  and  nonlinear  functions  are 
included  in  nonlinear  estimator  that  acts  on  the  model  regressors  to  provide  the  model  output  [9].  Figure  9  shows  the 
nonlinear  ARX  model. 


Nonlinear  estimator 


Figure  9:  NARX  Model  Structure 

The  algorithmic  steps  of  NARX  modelling  are  as  follows 

Step  1:  By  using  past  output  data,  present  and  past  input  values  calculate  the  regressor  value. 

Step  2:  By  using  the  nonlinearity  estimator  map  the  regressors  to  the  model  output. 

4.2.  Hammerstein  Modelling 

Hammerstein-Wiener  models  uses  linear  block  in  series  with  static  nonlinear  block  to  describe  the  dynamic 
systems  [10],  [11]  and  [12].  Figure  10  shows  basic  model  of  a  Hammerstein-Wiener  model. 
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Figure  10:  Model  Structure  of  Hammerstein-Wiener 

The  algorithmic  steps  used  in  this  modelling  are  as  follows 

Step  1:  By  using  input  data  calculate  w(t)  =f(u(tj) 

Step  2:  Using  initial  conditions  x(t)  =  (B/F)*w(t)  and  w(t)  calculate  the  output  of  the  linear  block  and  configure  it. 

Step  3:  By  transforming  the  response  of  x(t)  calculate  model  output  using  nonlinear  function  h  as  y(t)  =  h(x(t)) 
and  configure  it. 

In  most  of  the  estimation  methods  Prediction  Error  Minimization  (PEM)  method  is  used  for  error  minimization 
purpose.  The  algorithmic  steps  of  PEM  are  as  follows 

Stepl:  Choose  a  model  structure  of  the  form 

y{t)  =  B(p-l-0)*u{t)  +  A(p-x-0)*e{t),  E[  e(t)e(t)*  ]  =  MO)  (22) 

Step2:  Choose  a  predictor  of  the  form 

y(t  1 1-10) =Up~\O)*y(t)+M(p-\0)*u(t)  (23) 

Step3:  Select  a  criterion  function  / (R(0)) 

Step4:  Determine  6  that  minimizes  the  loss  function/ 

5.  MATHEMATICAL  MODELLING  OF  VRM 

The  difficulty  to  observe  internal  phenomena  and  the  complexity  of  the  VRM  system  do  not  allow  a  thorough 
knowledge  process.  Number  of  modelling  methods  are  attempted  to  overcome  these  difficulties  have  been  carried  out  in 
the  last  decades.  The  prediction  of  some  key  parameters  in  the  operational  process  of  VRM  is  very  important  for  its  safety 
and  reliability.  Because  it  works  too  long  in  poor  conditions,  its  key  operating  parameters  are  nonlinear  and  time-varying 
and  the  traditional  prediction  methods  are  difficult  to  achieve  high  accuracy  [13]  and  [14].  A  different  process  parameters 
available  VRM  process  is  inlet  feed,  inlet  pressure,  outlet  pressure,  mill  inlet  temperature,  mill  outlet  temperature,  hot  air 
supply,  water  injection,  Separator  speed,  mill  vibrations,  grinding  table  speed,  bed  thickness  and  power  consumption  of  the 
mill.  The  real  time  VRM  process  consists  of  three  main  and  three  supporting  rollers.  The  data  collected  VRM  mill 
specifications  are: 

(i)  OPC:  210tph,  26.5Kwh/MT  cement 

(ii)  PPC:  300tph,  19.8  Kwh/MT  cement 

(iii)  PSC:  170tph,  32.0  Kwh/MT  cement 

(iv)  Cement  grinding  specific  power:  24.5Kwh 
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(v)  Cement  mill  MTBF:  125Hrs 

(vi)  Fly  ash  absorption:  34.2% 

(vii) Slag  absorption:  58% 

From  the  industry  expert  knowledge,  the  selected  input  parameters  are  input  feed  (TPH),  water  injection  pump 
speed  (RPM)  and  output  parameters  are  mill  differential  pressure  (mbar),  mill  vibrations  (mm/sec),  mill  bed  thickness 
(mm)  and  mill  power  consumption  (Kwh)  for  system  identification.  Figure  11  shows  the  real  time  plant  data  used  for 
system  identification  after  pre-processing. 

6.  RESULTS 


From  Table  1  the  best  fit  model  for  the  process  variables  differential  pressure  and  mill  vibrations  is  given  by 
discrete  time  state  space  model;  for  the  process  variables  mill  bed  thickness  and  mill  power  consumption  the  best  fit  is 
given  by  the  method  ARMAX.  The  best  fit  models  are  shown  in  Table  2. 
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Figure  11:  Real  Time  Pre-Processed  Data  For  System  Identification 
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Table  1:  Comparison  of  Yarious  System  Identification  Methods  Applied  for  YRM  Modelling 


Modelliitg 

Metliod 

Model  Type 

Diffeientiiil 

Pressiue 

(inbar) 

°o  Fit 

Mill 

Yibrations 

(mmSec) 

°/o  Fit 

Mean 

Square 

Eiror 

Mill  Bed 

Tliickuess 

(lmn) 

0  o  Fit 

Mill  Power 
Consmnption 
(Kwh) 

%  Fit 

Mean 

Square 

Eiror 

State  space  model 

Continuous  time 

92.02 

89.1 

0.00135 

86.95 

91.67 

0.00291 

Discretetime 

92.58 

89.51 

0.001213 

87.29 

89.79 

0.003178 

Polynomialmodels 

LARX 

89.72 

88.04 

0.001906 

87.03 

92.36 

0.002761 

IV 

46.91 

35.5 

0.1413 

78.4 

79.6 

0.01044 

ARMAX 

89.44 

88.64 

0.001873 

88.53 

92.52 

0.002272 

OE 

51.98 

6.546 

0.07648 

28.91 

50.33 

0.09075 

BJ 

89.83 

81.27 

0.003174 

87.59 

92.14 

0.002621 

Non-  linear  models 

NARX 

89.61 

87.94 

0.001948 

88.14 

92.45 

0.002399 

Hammerstein-Wiener 

21.73 

21.84 

0.09685 

24.27 

72.54 

0.0806 

Transter  functionmodels 

Continuoustime 

45.13 

60.34 

0.05375 

60.81 

71.59 

0.02821 

Discretetime 

32.45 

42.16 

0.06465 

30 

24.48 

0.1243 

Table  2:  Best  Fit  Model  Results  for  VRM  Process 


Best  Fit  Modei  for  Differendal  Pressure  and  Mill  Yibrations  Using 
Discrete  State  Space  Model 

x(t  +  Ts  )  =  A  *  x(t)  +  3*  u  (t)  +  K  *  e(t ) 
y(t)  =  C*x(t)+  D*u(t)  +  e(t) 

Best  Model  for  MiD  Bed  Thickness  and  Mill 
Power  Consumption 

xi  x2  x*  x4  x5  x6  ar?  x%  x9 

»\  1016  -01776  -  195  J  -005647  12*7  -0022*  -06*67  -  07561  -02345 

*2  —  0*705  904*  -  .2112  1093  -  2231  07906  -  03992  3402  1463 

»3  15*1  06603  9174  .1149  -  01026  0*431  09*58  049*5  -  1024 

»4  1162  -  0471  -  04293  .4177  -  001098  2011  -  4061  -  3144  3033 

4".5  003797  2175  -  1378  3104  3CKT  4345  418'  1574  -1323 

*6  06148  -02107  -  0001*37  -  606  -  -267  3717  06327  -  5504  -  2526 

»7  -  1492  1135  09442  -  4456  -  2038  0141*  1193  012*9  3606 

x%  09623  01412  1*94  ,257*  1467  4261  -1  121  1761  -  273* 

>9  236  -  1171  0*264  -  5776  -.1237  -  J24J  *J06  0976  8397 

»2  -  06703  -  03855 

jt3  -  07161  .01397 

v  6  -  5421  1669 

m  7  -  ©  2974  O  3423 

«8  0  3149  0  01864 

c  -  .Vl  -  2  er  764-»  0-225  -  1966  013?  09712  0239'  1154  -  09566 

,v2  -  1  59«  -  1  652  08096  01398  05982  1426  1149  <0467  -  0766 

«#1  «2 

Z>  -  „vl  0  0 

yx  o  o 

'  3  4825  4512 

-  .  »4  -  .1102  -  2297 

*  5  -  243  -  1054 

Modtl  output  for  Mill  bcd  thidcwss 

4d*.v_«r)  - 

4(z>- 1-1. 406  r-*+ 0.4441  r*3, 

0.01827  0.01683  :'2 

B\(z)  =  -0.008054  :‘*+ 0.002487  . 

32(:)- 0.01578  :‘l +0.007271  :’J. 

C(:)-l  +  0.001888  r”  +  0 .9862  : "* 

Modd  output  for  Mill  powrr  consumption 
A(:)y_2(t)~-A_K:)y_Kt)^B(:)ii(t)^C(:)eJ(t) 

A(:)~  1- 1.646:'*  + 0.66461  :’2. 

.I_l(:)»  0.02539  :'1- 0.03069  :'2 
il(:)=0  004848  :'1  -0.01838  :*2. 

32(:)-  0.0009788  :'1  -  0.007837  :'2. 

C(:)-  1  +  0.09974  :"  + 0.1167  :" 

Thc  ordcrs  of  the  pohnomial  are:  na»[2  2}  2],  nb-[2  2;2 

2].  nc-[2;2],  nlc-[l  1;1 1] 

7.  CONCLUSIONS  AND  FUTURE  SCOPE 

Simulation  results  show  satisfactory  prediction  capabilities  for  VRM  modelling.  After  comparing  the  results  state 
space  and  Nonlinear  ARX  models  gives  better  fitted  models  when  compared  with  their  counterparts.  This  work  can  be 
extended  in  future  by  applying  methods  like  neural  network  time  series,  fuzzy  logic,  spectral  and  correlation. 
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